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We show that the onset of the snake instabihty of ring dark solitons requires a broken symmetry. 
We also elucidate explicitly the connection between imaginary Bogoliubov modes and the snake 
instability, predicting the number of vortex-anti-vortex pairs produced. In addition, we propose a 
simple model to give a physical motivation as to why the snake instability takes place. Finally, we 
show that tight confinement in a toroidal potential actually enhances soliton decay due to inhibition 
of soliton motion. 
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The nonlinearity of the equations of motion governing 
the behaviour of optical systems and ultracold atomic 
gases allows for interesting solutions such as bright and 
dark solitons [IHll , which remarkably propagate without 
dispersion [S] and scatter elastically. In nonlinear optics, 
this means that soliton light pulses sent through opti- 
cal fibers propagate without changing shape, whereas in 
Bose-Einstein condensates and superfluid Fermi gases [B] , 
matter wave solitons correspond to shape-maintaining 
dips or humps in the atomic density. 

In both settings dark solitons have been observed ex- 
perimentally [71 [5], but strictly speaking the stability 
holds only in one-dimensional systems. In higher dimen- 
sions, in general, they collapse into vortex- ant i- vortex 
pairs in 2D (or vortex rings in 3D) through the snake in- 
stability [To] , even though dark solitons retain many of 
their solitonic properties in nearly-integrable systems, for 
example, in the presence of an external trap |llj . The sta- 
bility and dynamics of dark solitons has been discussed 
theoretically 15j, in particular, complex frequencies 
in the Bogoliubov-de-Gennes (BdG) excitation spectrum 
have been demonstrated to drive the instability [IB] . 

Cylindrically symmetric systems, bringing forward the 
concepts of ring bright [17] and dark [T8H20] solitons 
(RDS), offer an example of a two-dimensional system, 
where the dynamics can be reduced to a one-dimensional 
equation, and thereby making the question of stability 
relevant. It has been observed that filling the notch of a 
RDS stabilises against the snake instability [5T] , a result 
we explain in this Letter. Specifically, ring traps [22H25] 
might stabilise the ring dark soliton, but we show that 
their effect is in fact the opposite. 

In this Letter we consider the stability of the RDS 
in cylindrically symmetric confinement, showing that it 
is stable as long as the symmetry in the ^-direction is 
maintained. We show that the type of the simulation grid 
plays an important role related to the induced breaking of 
this symmetry. We propose a physical model explaining 
why the snake instability takes place, and further eluci- 
date the connection between the complex BdG spectrum 
and the snake instability. Based on this model, we show 



how the number of vortex-anti-vortex pairs is directly 
related to the BdG spectrum. 

To investigate the properties of the two-dimensional 
condensate, we assume it is described by a macroscopic 
wavefunction, ip, which solves the Gross-Pitaevskii equa- 
tion in the (x, y) plane: 

iV;j, = -VV + KrapV' + C2I3|^|V, (1) 

where C2D = 4\/7rAfa/(aosc), where TV, a, and aosc are 
the number of atoms in the cloud, the s-wave scattering 
length of the atoms, and the characteristic trap length 
in the z-direction respectively. We have obtained dimen- 
sionless quantities by measuring time, length and energy 
in terms of oj~^ , a^sc = \ffij(^^™^x) and hio^ respectively, 
where l^^ is the angular frequency of the trap in the x- 
direction. This is equivalent to setting uj^ — h ~ 2to = 1. 
Similarly in the polar coordinates (r, 9) : 

i4>T = -0rr-^009+ ^^trap " 'Z'+^^I'/'P^, (2) 

where we have used the standard scaling (p — ^/rip. The 
normalisation of ip and (j) is to unity. Exact soliton-like 
solutions for Eq. ([2| are difficult to find, although some 
can exist for very specific forms of Vtrap [26) . 

The potential term Vtrap is the external trap, and we 
consider cylindrically symmetric traps, for which ujy — 
LOx = ijJr- In particular, we consider harmonic ring traps 
with toroidal potentials 

V,,,^ = -^{R~Rs)\ (3) 

where Rs is the radius of the trap minimum, and fc is a 
real nonzero constant. 

To find the ground and excited states and to consider 
time evolution, we propagate Eqs. ([T]) and ^ in imag- 
inary and real time respectively, using a split-operator 
Fourier method [57]. We consider Cartesian {x,y) and 
polar (r, 9) coordinate meshes and for the polar mesh we 
extend the grid to the negative r-axis and demand the 
wavefunction be antisymmetric about r = 0. The ring 







solitons are generated by printing an appropriate phase 
step in the radial direction. 

We also explore the Bogoliubov quasiparticle (or 
normal eigenmode) excitation spectrum, defining the 
Bogoliubov amplitudes u and v by {u,v}{r) — 
e^'^^{uq,Vq}{r) [55] representing a partial wave of an- 
gular momentum q relative to the condensate. The 
Bogoliubov-de-Gennes equations for low-lying modes Q,q 
in dimensionless harmonic oscillator units {uj^ — h — 
2m = 1) are 



n c. 



n 



(4) 



where U = -V^ g -I- Vtrap + 2C2_d|(/>oP - and 0o is 
the cylindrically symmetric self-consistent wavefunction, 
containing a ring dark soliton. Here /i is the chemi- 
cal potential, and u and v satisfy the general condition 



Vq) = 0. The instability time 



(f2g — Qqi) J dr[UqUqi — uqf uq 

is given by 1/Im(r2q), and if all eigenvalues are real, the 
condensate is dynamically stable. We employ Lagrange 
functions [29, together with Hessenberg QR iteration [30] . 

We have observed the snake instability of the ring dark 
soliton only when the symmetry in the ^-direction is bro- 
ken. This can happen, for example, by using a rectangu- 
lar coordinate system (x, y) and a symmetrical potential, 
or artificially adding an anisotropy in the trapping poten- 
tial and using a (symmetrical) polar coordinate system 
(r,9). In particular, in our simulations the ring dark 
soliton is a stable object as long as this 0-symmetry is 
maintained. The previous results with clear decay into 
a vortex-anti-vortex necklace in the literature [19] follow 
because a rectangular grid was used. It is clear that in 
such a situation the vortex-anti-vortex pairs also appear 
in multiples of four because the simulation setup has re- 
flection symmetry about two orthogonal axes (what hap- 
pens in one quadrant is copied 4- fold). 

When there exists symmetry in the ^-direction, the 
Gross-Pitaevskii equation governing the behaviour of the 
system becomes effectively one-dimensional in the radial 
direction. In the rectangular case, for the equation to 
become effectively one-dimensional, one would need an 
infinite system that is homogeneous in one direction (or 
in 3D homogeneous in x and y if the dark soliton is in 
the z-direction) . In Ref. [T^] it was shown that such a 
planar dark soliton is stable as long as the transverse di- 
mensions are restricted enough to make the system quasi- 
one-dimensional, but a similar quasi-one-dimensionality 
is achieved by the transversal symmetry, i.e. in the zero 
limit of the trapping. When the symmetry is broken, the 
snake instability starts to develop at the location(s) of the 
symmetry breaking. Note that the symmetry breaking is 
not spontaneous, but induced. 

On the torus geometry, we have in general observed 
the snake instability provided we use polar coordinates 
and an asymmetric potential (see Fig. [ij or rectangular 
coordinates (see Fig. [2|, but the trap parameter k, the 
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FIG. 1. Snapshots of real-time propagation of the GPE ([2| 
in a polar grid at (a) t = 0.0, (b) t — 1.5, (c) t = 3.0, and 
(d) t = 4:A/ojx with the harmonic torus ([3| with k = 96, 
Rs ~ 8.5, and C2D ~ 400, but with a slight added anisotropy 
a,t 9 = n/2. The transverse oscillations of the snake insta- 
bility are seen to start immediately originating at the trap 
anisotropy and propagate along the ring with a finite speed, 
in accordance with the model discussed in the text. 
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FIG. 2. Snapshots of real-time propagation of the GPE ([T]) 
in a rectangular grid at (a) t = 0.0 and (b) t — with 
the same parameters of Fig. [l] excluding the anisotropy. The 
transverse oscillations of the snake instability are seen to orig- 
inate over the whole ring resulting in the collapse of the ring 
dark soliton (producing eventually a necklace of vortex-anti- 
vortex pairs). 



nonlinear coupling constant 6*213, the soliton radius Rs, 
and the particular way in which the symmetry is broken 
affect the decay time. To eliminate motional dynamics we 
take into account the effective curvature-induced poten- 
tial, and provided the radius of the ring is large enough, 
the stationary radius is determined by |19j 



^(^trap)/? + ^ 



0. 



(5) 
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Our simulations show that the higher the curvature of the 
trap (in the radial direction), the faster is the onset of 
the snake instability. Also, the decay time, T^, decreases 
if we increase the nonlinearity while keeping the other 
parameters fixed. In a fixed harmonic trap ([s]) with k = 1 
and Rs = 14.5, we find that Td - C:^^^^. 

For a fixed nonlinearity of C2D = 400, the ring dark 
soliton survives for « 10 times longer in a harmonic trap 
(i?5 = 0) than in a torus of similar size. The survival 
times for the Rs — case are much longer because the 
ring dark soliton can move quite freely, whereas toroidal 
confinement restricts motion, which leads to faster de- 
cay. If the soliton moves fast enough, it will fade out 
without the snaking instability. The effect of the motion 
is in line with the results of the experiments for moving 
planar dark solitons [3 [3]. Furthermore, if the soliton 
is stationary, also the vortex-anti- vortex necklace will be 
stationary. 

We have also observed periodic revival and subsequent 
decay into a vortex-anti-vortex necklace of the ring dark 
soliton in the harmonic torus ([3| with k = 1, Rs — 7.5, 
and C2D = 400. The initial radius of the soliton is given 
by Eq. The initial ring dark soliton decays to a 

vortex-anti-vortex necklace, but is revived several times 
later on with a period of « 8.8 /uj^- 

One should note that tightening the trap or increas- 
ing the nonlinearity have a similar effect on the healing 
length. 
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where no is the bulk density. It characterises both the 
width of the density notch in a dark soliton as well as the 
size of vortex cores. Also, ^ will decrease if we make the 
trap tighter (to preserve normalisation) or increase the 
nonlinearity. With decreasing healing length it becomes 
easier for the snake instability to deform the condensate 
with respect to the bulk. 

Of course, nothing would happen if there was no driv- 
ing force trying to deform the condensate around the soli- 
ton's notch in the first place. As a first approximation, 
this force can be understood if we replace a homogeneous 
bulk condensate with an infinite potential for x < xs, 
where xs marks the position of the wall. The infinite 
wall forces a zero boundary condition for the wavefunc- 
tion just like a planar dark soliton at 2:5, so for x > xs 
nothing changes. This happens because the kink 



ip — tanh 



(x > 0) 



(7) 



solves Eq. ([TJ {C2D = 1) with an infinite wall potential 
at a; < and a zero potential elsewhere (uniform for y 
and z). Therefore, the dark soliton strip can be thought 
of as serving the role of a two-sided wall, one side facing 
towards increasing coordinate and the other side in the 




FIG. 3. The hydrodynamical force (solid), Eq. Kt, around the 
notch of a planar dark soliton at 2; = (dashed). Notice the 
sign of the force and that the force scales as ~ explaining 
the relationship between the decay time and ^. 



opposite direction. Sudden removal of such a wall poten- 
tial will result in the expansion of the condensate to the 
previously forbidden area. This can be seen from the hy- 
drodynamical description of the condensate {tp — \fne^' 
where S is the phase): 
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dt 
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(8) 

where -p is the pressure, v = ^ V5 the fluid velocity, and 
n is the (local) density. Since the length scale of spatial 
variation of the condensate wavefunction at the soliton's 
notch is of the order of ^, the relevant and dominant term 
is the quantum pressure. Focusing on the planar case, 
and measuring length, time, and density in terms of ^, 
and no /TV, respectively = h = 1m = uqIN = \) 
the quantum pressure term in Eq. (|8|) reduces to 



1 dv 
2 



1 



(9) 



where the extra minus sign follows from the soliton wave- 
function {S — > —5). Equation ([9]) is plotted for a dark 
soliton strip n = tanh^ [x j \/2) in Fig. [s] The force is 
pointing towards the notch on both sides, and in partic- 
ular, it scales as ~ (T^ ■ Filling the notch with a repulsive 
component reduces the effect of this force, explaining 
why it acts as a stabilising factor. 

Provided this force can point at least infinitesimally 
askew, the two sides will then eventually slip past each 
other, twisting the density notch of the soliton as they 
expand. Since there is a phase difference of tt between 
the two sides, twisting the density notch by such an ex- 
pansion will result in an accumulation of phase by 27r 
when we go around one dot-like density notch the twist- 
ing leaves, i.e. vortices. The phase singularity generates 
a topological protection for the perseverance of the vor- 
tices. This process will then propagate along the ring 
(see Figs. [1] and [2|. 

The number of vortex-anti-vortex pairs is then deter- 
mined by the length of the soliton and the size of the 
vortex core. As a first approximation we can say that 
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FIG. 4. The correspondence between the predicted (dashed, 
see Eq. (lOl) and measured (solid) number of vortex-anti- 
vortex pairs for the harmonic potential ([3| with k — 1 and 
C2D = 400. The healing length was measured by finding a 



best fit to 



/ R-Rs \ 



For the low radii the 



ring soliton disappears as sound waves without vortices. The 
shaded areas show where the complex Bogoliubov modes such 
that g is a multiple of 4 are active. In the middle of the regions 
the modes are also the fastest ones (see Fig. [5|, but they 
are still chosen by the condensate near the region boundaries 
where other modes would be faster as such, because A''„ must 
be a multiple of four by symmetry reasons. 



the number of vortex-anti- vortex pairs, Ny, will be 



2ttRs 



TT^aUQ 
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(10) 



naturally such that Ny is an integer (see Fig. |4]) because 
each vortex-anti- vortex pair and the distance to the next 
pair needs a total length of ~ 8^ on the circumference. 
Their locations within the ring is determined by where 
the symmetry was broken. 

The Bogoliubov decay times define a quantitative time 
scale related to the snake instability (see Fig.js]). Accord- 
ing to these times, a larger radius leads to a longer decay 
time, which can be understood because in a larger trap 
the bulk density is reduced to maintain normalisation. 
This increases ^ and hence the decay time. 

Because we have observed the snaking instability in 
general, and because there is only one imaginary Bogoli- 
ubov mode (per g), it is reasonable to assume they are 
related. The connection between the form of the com- 
plex Bogoliubov modes and the planar snake instability 
has been discussed in Ref . |16| , and here we report similar 
results for the ring and elucidate further the role of the 
complex modes. Firstly, we note that for an exact soliton- 
like dark ring solution of Eq. ([!]) , the complex Bogoliubov 
modes are very similar to Fig. [5] and show transverse 
oscillations around the ring resembling the snake insta- 
bility 26^ . This is the case for the numerically imprinted 
ring dark soliton as well (see Fig. [6]). 

In particular, the number of vortex-anti-vortex pairs, 
Ny, is determined by the decay channels, the channel q 
giving q pairs (see Fig. |4]) . Due to the symmetries of the 
rectangular square grid, Ny must be a multiple of four, 
but also the fastest decay channels (i.e. largest imaginary 



eigenmodes) will be preferred. Furthermore, Eq. ( 10 ) sets 
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FIG. 5. Results of the Bogoliubov analysis for the harmonic 
potential ([3| with fc = 1 and C2D = 400 such that the soli- 
ton is printed at the radius given by Eq. ([5|. Shown are the 
Bogoliubov eigenvalues for modes with a dynamical insta- 
bility, labelled by the value of q. There is only at most one 
imaginary eigenvalue per q. The inset shows the reciprocal 
numbers giving the decay times. The solid black line with- 
out markers is given hy Q = 1.0/(0. 135i? + 0.444). For low 
radii the primary decay channels of the RDS are the q = 2, 3 
modes, while higher modes are dynamically stable. As the 
radius is increased the higher modes become unstable as well. 
The shading is the same as in Fig.|4] When there is no soliton 
all the eigenvalues are real. 




an (approximate) upper bound, and of course the way in 



FIG. 6. (Colour online.) The wave functions |</>oP (solid) and 
\4>o + u + (dashed) at 6 = 2tv/8 for a ring dark soliton in 
the harmonic trap ([3| with fc = 1 and C2D = 400 such that 
the soliton is printed at _R = 8.65407 (as given by Eq. ([5|). 
u and V correspond to the only imaginary Bogoliubov eigen- 
value — 0.613i of the mode q = 8. The inset shows how the 
location of the soliton's notch is affected by the Bogoliubov 
mode around the circumference, mimicing the snake instabil- 
ity and coinciding exactly with the locations of the vortices 
(the Bogoliubov mode in the inset is exaggerated for clar- 
ity). For this radius, the g = 8 mode is the fastest, and it 
produces q = Ny = 8 pairs (see Fig. [4|, in accordance with 
observations. 



which the symmetry is broken dictates how the snake 
instability can take place. The exact choice of qs is then 
a balance of several factors. 

In summary, we have studied the stability of ring dark 
solitons on the torus geometry, elucidating the connec- 
tion between relevant symmetries, the BdG spectrum and 
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the snake instability. In addition to assessing the role of 
toroidal confinement in ring soliton decay times, we have 
provided a physical description for the decay of dark soli- 
tons into vortices. Even though we have considered a ring 
dark soliton, it is important to note that our modelling 
is general in nature. The location(s) of the symmetry 
breaking are seen to be the origins of the snake instabil- 
ity in many other configurations as well [9l 1311 132j . 
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